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A simple test for periodic signals in red noise 
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Abstract. We demonstrate a simple method for testing the significance of peaks in the periodogram of red noise data. The 
procedure was designed to test for spurious periodicities in X-ray light curves of active galaxies, but can be used quite generally 
to test for periodic components against a background noise spectrum assumed to have a power law shape. The method provides 
a simple and fast test of the significance of candidate periodic signals in short, well-sampled time series such as those obtained 
from XMM-Newton observations of Seyfert galaxies, without the need for Monte Carlo simulations. A full account is made of 
the number of trials and the uncertainties inherent to the model fitting. Ignoring these subtle effects can lead to substantially 
overestimated significances. These difficulties motivate us to demand high standards of detection (minimum > 99.9 per cent 
confidence) for periodicities in sources that normally show red noise spectra. The method also provides a simple means to esti- 
mate the power spectral index, which may be an interesting parameter itself, regardless of the presence/absence of periodicities. 
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1. Introduction 

Many astrophysical sources show erratic, aperiodic brightness 
fluctuations with steep power spectra. This type of variability is 
known as red noise. By 'noise' we mean to say that the intrinsic 
variations in the source brightness are random (this has nothing 
to do with measurement eiTors, also called noise). Examples 
include the X-ray variability of X-ray binaries (XRBs; e.g. van 
der Klis' 1995t and Seyfert galaxies (e.g. Lawrence et al. 1987 
Markowitz et al. [2003 ). The power spectrum of these varia- 
tions, which describes the dependence of the variability ampli- 
tude on temporal frequency, is often reasonably approximated 
as a simple power law (over at least a decade of frequency). 
This featureless continuum spectrum does not offer any char- 
acteristic frequencies that could be used as diagnostics. 

XRBs often show quasi-periodic oscillations (QPOs) that 
show-up as peaks in the power spectrum over the continuum 
noise spectrum. These can be thought of as half-way between 
strictly periodic variations (all power concentrated at one fre- 
quency) and broad-band noise (power spread over a very broad 
range of frequencies). A combination of periodic oscillations 
with similar frequencies, or a single oscillation that is per- 
turbed in frequency, amplitude or phase can produce a QPO. 
QPOs are one of the most powerful diagnostics of XRB physics 
(see e.g. van der KHs 1995. M'Clintock & Remillard 2004 1. 
The detection of periodic or quasi-periodic variations from 
a Seyfert galaxy would a be a key observational discovery, 
and could lead to a breakthrough in our understanding if the 
characteristic (peak) frequency could be identified with some 



physically meaningful frequency. For example, if we assume 
a 1 /Mbh scaling of frequencies we might expect to see ana- 
logues of the high frequency QPOs of XRBs in the range 
/qpo ~ 3 X 10-^(Mbh/10^Ms)-^ Hz (Abramowicz et aL l2004t . 

However, claims of periodic variations and QPOs in the X- 
ray emission of Seyfert galaxies have a chequered history, with 
no single example withstanding the test of repeated analyses 
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and observations (see discussion in Benlloch et al. l2001t . The 
confusion arises partly due to the lack of a standard technique 
to assess the significance of a periodicity claim against a back- 
ground assumption of random, red noise variability. Indeed, 
as Press (1978) and others have remarked, there is a tendency 
for the eye to identify spurious, low frequency periods in ran- 
dom time series. Tests for the presence of periodic variations 
against a background of white (flat spectrum) noise are well 
established, from Schuster ( 1898) and Fisher J 19291 1. these are 
reviewed in section 6.1.4 of Priestley il98U . and discussed in 
an astrophysical context by Leahy et al. (|1983) and van der 
Klis (1989 1. But without modification these methods cannot 
be used to test against red noise variations. Timmer & Konig 
( 1995) and Benlloch et al. J200H have proposed Monte Carlo 
testing methods applicable to red noise but the relatively high 
computational demands of these methods may be enough to 
deter some potential users. Israel & Stella ( 1996 1 proposed a 
method that does not require Monte Carlo simulations but is 
not optimised for short observations of power law continuum 
spectra. 

This paper puts forward a simple test that can be used to 
test the significance of candidate periodicities superposed on 
a red noise spectrum which has an approximately power law 
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shape. The price of simpHcity, in this case, is that the test 
is only strictly valid when the underlying continuum spec- 
trum is a power law. The basic steps of the method are: (/) 
measure the periodogram, (//) estimate the red noise contin- 
uum spectrum and (///) estimate the significance of any peaks 
above the continuum. The stages of the method are explained 
in detail in the following sections. Section |2] gives a brief 
introduction to the statistical properties of the periodogram. 
Section |3] discusses a simple method for estimating the pa- 
rameters of a power law-like spectrum and sectionl^discusses 
how to estimate the significance of a peak above the contin- 
uum. Section |5l then demonstrates the veracity of the method 
using Monte Carlo simulations and section|6lreviews some im- 
portant caveats that must be considered when using this (and 
other) period-searching methods. Finally, section0gives a brief 
review of the method in the context of observations of active 
galaxies. The appendix discusses a more generally applicable 
method of periodogram fitting (that makes no assumption on 
the form of the underlying spectrum). 

2. The periodogram 

Given an evenly sampled time series Xk of K points sampled at 
intervals AT" we can measure its periodogram (Schuster 189^, 
which is simply the modulus-squared of the discrete Fourier 
transform, X(fj), at each of the n - K/2 Fourier frequencies: 

2AT , 

- ^1^/ (1) 

The normalisation is chosen such that the units of the peri- 
odogram are (rms/mean)^ Hz"' (where rms/mean is dimension- 
less) and summing the periodogram over positive frequencies 
gives the sample variance in fractional units. The periodogram 
is evaluated at the Fourier frequencies fj = j/KAT with 
j - 1,2, . . . , K/2. The original purpose of the periodogram 
was as a tool for identifying 'hidden periodicities' in time se- 
ries. However, the periodogram of a noise process, if measured 
from a single time series, shows a great deal of scatter around 
the underlying power spectrum. Specifically, the periodogram 
at a given frequency, /(/,), is scattered around the true power 
spectrum, 'Pifj), following a distribution with two degrees 
of freedom: 

/(/,) = nmla, (2) 

where ^2 is a random variable distributed as with two de- 
grees of freedom, i.e. an exponential probability distribution 
with a mean and variance of two and four, respectively: 

p^2(x) = e-''/2/2 (3) 

The periodogram is distributed in this way because the real 
and imaginary parts of the DFT are normally distributed for a 
stochastic process', and the sum of two squared normally dis- 
tributed variables is a ;if2 "distributed variable (Jenkins & Watts 
[T9681 Priestley [T98T1 Chatfield 1989 Bloomfield 2000 1. See 

' The DFT at the Nyquist frequency is always real when N is even 
so the periodogram at this frequency is distributed asx'[, i-e. with one 
degree of freedom. 



also Scargl e iTWh . Leahy et al. ( fT^ . Wall & Jenkins (I^UU^ 
and Groth ( I1975> for further discussion of this point. 

If the spectrum is flat ('white noise') and its power level 
known a priori then we can simply make use of the known 
probability distribution (equations|2and|3} to estimate the Uke- 
lihood that a given periodogram ordinate exceeds some thresh- 
old. If the power level is not known there is an added un- 
certainty, but nevertheless an exact test does exist (Fisher's g 
statistic: Fisher 1929; section 6.1.4 of Priestley 1981 ) to esti- 
mate the likelihood that the highest peak in the periodogram 
was caused by a random fluctuation in the noise spectrum (see 
alsoKoen 1990). 

For the more general case of non-white noise there is no 
such exact test. When examining the periodogram of red noise 
data such as from Seyfert galaxies, we need to be careful not 
to identify spurious peaks. Even in the 'null' case (i.e. no pe- 
riodic component) peaks may occur in the periodogram due to 
sampling fluctuations. In particular the eye may be drawn to 
low frequency peaks because, in red noise data, there is much 
more power and more scatter in the periodogram at low fre- 
quencies. Given a large amount of data we can average the pe- 
riodogram in one of the standard ways (see e.g. van der Klis 
0.989; Papadakis & Lawrence 1 1993a> . fit the continuum using 
a standard ;t'^-minimisation tool (e.g. Bevington & Robinson 
IT9921 Press et al. lT996.) and test of the presence of addition fea- 
tures. This is the standard procedure for analysing XRB data. 
If we have a very limited amount of data, such that we cannot 
afford to average the periodogram, we are faced with a more 
difficult situation. Figure^gives an example like this. The pe- 
riodogram of a short time series, containing red noise (gener- 
ated using the method of Timmer & Konig ,1995J . shows a large 
peak at / = 4 X 10"^. Could this be due to a real periodic varia- 
tion present in the data or is it just a fluctuation in the red noise 
spectrum? 

3. Fitting the periodogram 

3.1. Least squares (LS) fit to log-periodogram 

If the underlying power spectrum is suspected to be a power 
law then then parameters of interest are its slope, a, and nor- 
malisation A^. One of the simplest methods to estimate these pa- 
rameters from the raw (unbinned) periodogram is to fit it with 
a model of the form !P(/) - Nf^" using the method of least 
squares (LS). The problem with this is that the periodogram 
is distributed around the true underlying spectrum in a non- 
Gaussian fashion and, more seriously, the distribution depends 
on the spectrum itself (equation|5J. 

To simplify the problem we can fit the logarithm of the pe- 
riodogram, as discussed in some detail by Geweke & Porter- 
Hudak (il983> see also Papadakis & Lawrence .1993a< . The 
scatter in the periodogram scales with the spectrum itself; the 
scatter is multiplicative in linear-space. This means the scatter 
is additive in log-space: 

log[/(/;)] = log[;P(/>)] + \og\xll2] (4) 

and therefore identical at each frequency (the data are ho- 
moskedastic). Working with the logarithm of the periodogram 
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Frequency f 

Fig. 1. Periodogram of a short (K = 256) time series contain- 
ing red noise. The upper panel shows the periodogram using 
linear axes, the lower panel shows the same data on logarith- 
mic axes. The periodogram shows a red noise spectrum rising 
at lower frequencies. But the periodogram also shows a peak at 
f - Ax 10"^. Is this due to a real harmonic (periodic) variation 
or an artifact of the fluctuating noise spectrum? 

has two further advantages. The first is that if the power spec- 
trum is a power law, the natural way to plot the data is in log 
space. The power law becomes a linear function: log[!P(/)] = 
log[A^] - Q'log[/]. The second advantage is that the distribution 
of periodogram ordinates becomes less skewed. This reduces 
the effect of 'outliers' on the fitting. 

We must be careful fitting the logarithm. The expectation 
value of the logarithm of the periodogram is not the expectation 
value of the logarithm of the spectrum. However, the bias is a 
constant (due to the shape of the ;^'2-distribution in log-space) 
that can be removed trivially: 



<log[/(/,)]) = (logTO)]) + {\og\x\l2]) 



(5) 



Using equation 26.4.36 of Abramowitz & Stegun (I1964t we 
have {\og\xll2]) = -0.57721466 .../ ln[10]. (The number in 
the numerator is Euler's constant.) Therefore 



<logTO)]> = <log[/(/,)]) + 0.25068 . 



(6) 



The logarithm of the periodogram ordinate, with the bias re- 
moved (i.e. the 0.25068 added^), is thus an unbiased estimator 
of the logarithm of the spectrum, and is distributed indepen- 
dently and identically (about the underlying spectrum) at each 
frequency. (The raw periodogram points are not identically dis- 
tributed in linear space since the scatter depends on the spec- 
trum, which is a function of frequency.) Thus we can use a 
LS fitting procedure to get a reasonable estimate of the power 
spectral slope a and normalisation by fitting a linear func- 
tion y - mx -1- c to the plot of log [/(/,■)] versus log[//]. The 
slope of the linear fit gives a = -m and the y-intercept gives 
log(A^) = c + 0.25068. These give our estimate for the contin- 
uum: Pj = NfJ", or equivalently \og[Pj] - log[A] -a\og[ fj]. 

^ This is a more precise approximation to (log|;^^5/2]> than used by 
Papadakis & Lawrence tl993a< . 




slope estimate a 




Norm, estimate log[N] 

Fig. 2. Distribution of the slope and normalisation estimators 
derived from 10^ Monte Carlo simulations of K - 256 point 
time series (histogram). The 'true' spectral parameters were 
a - 2 and log[A] = 0. The predictions of Gaussian uncer- 
tainties, with widths given by equations ™d |8] are shown 
with the smooth curves. 



It is important to note that the datum at the Nyquist fre- 
quency (j - n) should be ignored in the LS fitting. This is 
because, as mentioned previously, the distribution of the pe- 
riodogram ordinate at this frequency is not identical to that 
at other frequencies (it follows a x\ distribution). This mi- 
nor detail means that the LS fit should be performed on the 
«' - n-l lowest frequencies that are identically distributed (in 
log-space). 

A drawback of fitting the periodogram, rather than the 
binned or averaged periodogram, is that it does not provide a 
in-built goodness-of-fit test. By binning the periodogram (van 
der Klis 119891 Papadakis & Lawrence I1993a> we can obtain 
Gaussian errors on each ordinate to be used in a;t'^-test. We do 
not have Gaussian error bars for the unbinned log -periodogram. 
But, since we know the expected distribution of the peri- 
odogram ordinates about the true spectrum we can compare 
this to the distribution of residuals from the fitted data using 
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a Kolmogorov-Smirnov test (Press et al. '19961. Specifically, 
we can compare the data/model ratio (in linear space) given by 
jj - 2Ij/Pj with the theoretical x^ distribution, if the model 
is reasonable yj should be consistent with the distribution. 
Furthermore, the KS test is most sensitive around the median 
value, and less sensitive at the tails of the distribution, which 
means that even in the presence of a real periodic signals (i.e. 
a few outlying powers) the test should give a good idea of the 
overall quaUty of the continuum fit. 

3.2. Uncertainties on the parameters 

The uncertainties on the slope and normalisation estimates 
from the LS method can be derived using the standard theory 
of linear regression (e.g. Bevington & Robinson 1992 Press 



et al. I1996 >. The eiTor on the slope (index) and intercept (log 
normalisation) are: 



err [a] 



and 



err^[\og(N)] = 
where aj = log[/y] and 

// f n' 

>=i Vj=l 



(7) 



(8) 



(9) 



and also cr^ - 7r^/6(ln[10])^ is the variance of the log- 
periodogram ordinates about the true spectrum (Geweke & 
Porter-Hudak il983t . The covariance of the two parameter esti- 
mates is given by: 



cov{a, \og[N]) 



(10) 



Here n' is the number of frequencies used in the fitting. 
Normally n' - n - \ since only the Nyquist frequency is ig- 
nored (because the periodogram at the Nyquist frequency does 
not follow the same distribution as at the other frequencies). 

The accuracy of these equations was tested using a Monte 
Carlo simulations. An ensemble of random time series, each of 
length K, was generated (using the method of Timmer & Konig 
11995 '!. For each series the power spectral slope and normal- 
isation were estimated using the LS method discussed above. 
Figure|2shows the distribution of estimates for 10^ realisations 
of time series generated by a process with an a = 2,A^ = 1 
spectrum. With only n - 127 periodogram points (i.e. K - 
256) the distribution of the estimates is reasonably close to 
Gaussian. (The distribution of .^V is log-normal because the esti- 
mated quantity log[.;V] is normally distributed in the LS fitting.) 
These two parameters are covariant in the fit; a low estimate of 
the slope tends to be correlated with a high estimate of the nor- 
malisation. Figure|3]illustrates the covariance between the two 
estimated parameters. The shape of these distributions is inde- 
pendent of the spectral slope, this was confirmed using Monte 
Carlo simulations of spectra with slopes in the range a = - 3. 




-0.2 0.2 

Norm, estimate log[N] 

Fig. 3. Demonstration of the covariance in the estimates of 
slope and normalisation from from LS fitting to the logarithm 
of the periodogram. The plot shows the results of fitting 5, 000 
Monte Carlo simulations of an a = 2, = 1 spectrum. 

The uncertainties on a and log[A], and their covariance, 
were estimated for different length series by the same Monte 
Carlo method as discussed above. These Monte Carlo uncer- 
tainties compare well with the theoretically expected uncertain- 
ties for the LS fitting method as discussed above (Fig.0}. 



3.3. Uncertainty on the model 

The expected uncertainties and covariance of the two model 
parameters can be combined to give an estimate of the uncer- 
tainty of the logarithm of the model, \og[Pj\, at a frequency fj, 
using the standard eiTor propagation formula. 



err^{\og{P{fj)}] = err\a] X {\og{fj]f + err^ [log(A)] 
2cov[a,log(A^)]x(log[/,]) 



(11) 



The first term in the sum accounts for the uncertainty on the 
slope (equation Q, the second accounts for the uncertainty in 
the normalisation (equation|8} and the third accounts for their 
covariance (eauationllO>. The uncertainty is frequency depen- 
dent. Even in log-space the model is more uncertain at lower 
frequencies, this simply reflects the fact that there are many 
more points at high log[/] than at low log[/]. But also note 
that the uncertainty on the logarithm of the model is indepen- 
dent of the model itself (slope and normalisation). 

The distribution of the power in the model (in log-space), 
log[!Py], is expected to be Gaussian with a width determined by 
the formula above. In linear-space the uncertainty on the model 
power, Vj, is log-normally distributed. The probability density 
function for the model power is therefore 



Ppp) 



1 



exp < 



{\n[y]-Mjf 
252 



(12) 



where Mj - ln[!Pj] is the expected value of the power (in 
log-space) and S j is the rms width of the distribution of pow- 
ers (also in log-space) as given by equation ^2 Both of these 
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Fig. 4. Demonstration of the uncertainties from LS fitting to 
the logarithm of the periodogram. For each value of K, 10^ 
time series were simulated (with an or = 2, = 1 spectrum). 
The periodogram of each time series was fitted (in log space) 
with a linear model. From the 10^ estimates of the two param- 
eters (slope and normalisation) their rms and covariance were 
measured (solid squares). The solid lines mark the predictions 
of equationslTllSlandlTUl 

are frequency dependent. Note that the log-normal distribu- 
tion is conventionally defined in terms of the natural logarithm, 
whereas previously the results were given in terms of base 
10 logarithm. The uncertainty on the model log -powers from 
equation[nneeds to be corrected: 

Sj^err[\og{P{fj)]]x\n[\0] (13) 

Figure|5]shows the distribution of model powers (at two differ- 
ent frequencies) for 10^ Monte Carlo simulations of K - 256 
time series. Clearly the predicted log-normal distribution (with 
parameters Mj and S j as given above) gives a good description 
of the real uncertainty in the model. 

3.4. Bias in the parameters 

Although in general the LS method does not yield the max- 
imum likelihood solution for non-Gaussian data, for the spe- 
cific case of a power law spectrum the parameters obtained 
from the log-periodogram regression, namely a and log[.;V], 
are unbiased. Figure |6l demonstrates this using Monte Carlo 
simulations. However, it should be noted that because the pa- 
rameter log[.;V] is normally distributed the parameter N will be 
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Model power 

Fig. 5. Monte Carlo demonstration of the distribution of power 
in the model (power law) spectrum. Using an a = 2, A^ = 1 
spectrum, 10^ random time series of length K = 256 were gen- 
erated. For each one the periodogram was fitted as discussed 
in the text. The histograms show the distribution of the power 
in the resulting 10^ models at the /' = 10 and j = 100 Fourier 
frequencies (fj = j/KAT). The solid line shows the predicted 
log-normal distribution (eauationll2t. 

log-normally distributed. Thus the mean value of .^V is not a 
good estimator of the true value (it will be biased upwards due 
to the long tail of the log-normal distribution). 

3.5. Summary of LS fitting 

The following summarises the LS fitting method: 

- Measure the periodogram of the time series. 

- Ignore the Nyquist frequency point. 

- Take the logarithm. 

- Fit a straight line using the standard LS method. 

- Test the goodness-of-fit using the KS test. 

- Estimate uncertainties on the parameters and the model. 

The uncertainties and covariance of the parameters can be 
computed using equations 171 Isl and 1 1 01 These can then be used 
to calculate the uncertainty on the logarithm of model power 
S j (equationO. 

4. Confidence limits 

4.1. The ideal case 

If we know the exact form of the spectrum we can divide this 
out of the periodogram. From equation |2l we can see that the 
ratio -yj = 2I{fj)/'P{fj) will be distributed like;^'^. We can use 
our estimates a and .^V to form the null hypothesis: the data 
were generated by a process with a spectrum P{fj) = ■^Z/" 
and no periodic component. We can estimate the probability 
that a large peak will occur in the periodogram, assuming the 
model spectrum, by comparing fj to the xl PDF (Priestley 
1981;Scargle 1982). 



6 



S. Vaughan: A simple test for periodic signals in red noise 



o 
oi 



8 w 



o 
d 



o - 



o 
d 

I 



—1 1 — r- 



I 1 1 1 1 1 1 — I I I 1 1 1 1 1 1 — I I I 1 1 1 1 1 1 — \—\- 



10 



100 



N 



1000 



Fig. 6. Demonstration of bias in LS fitting of the logarithm of 
the periodogram. For each value of K, 10^ time series were sim- 
ulated (with an a - 2, N = 1 spectrum) and the periodogram of 
each time series was fitted (in log space; see text) with a linear 
model. The parameters estimates were averaged over the 10^ 
realisations. The estimators a and log[./V] are clearly unbiased. 

We can define a (1 - e)100 per cent confidence limit on y, 
call this y^, as the level for which, at a given frequency, the 
probability of obtaining a higher value by chance is Pr{y > 
7^) = e on the assumption that the null hypothesis is true. The 
chosen value of e represents the 'false alarm probability'. The 
integral of the xl probability distribution gives the probability 
of a single sample exceeding a value of by chance: 



Prix > Je 



1 r 

2 Jr. 



--V/2 



dx ■■ 



For a given probability e, we can rewrite this: 



-2 ln[e]. 



(14) 



(15) 



For example, using e = 0.05 (i.e. a 95 per cent significance test) 
we find 70 05 = 5.99. This means that if the null hypothesis is 
true the probability of the ratio jj being higher than 5.99 is only 
0.05. We can therefore define our 95 (and 99) per cent confi- 
dence limits on the log-periodogram as the model P{f) - Nf^" 
multiplied by the appropriate 7^/2. (In log-space we simply add 
the appropriate log[7f/2] to the model.) 

However, these confidence bounds correspond to a single 
trial test, they give the probability that a periodogram point at 
one particular frequency will exceed jeV Usually there are 
n' — n- \ independent trials since only the Nyquist frequency 



is ignored (leaving n' independently distributed periodogram 
points to be examined). We must account for the number of 
independent trials: 



7, = -21n[l-(l-e„0'^"'] 



(16) 



This gives the value of y that has a probability of being ex- 
ceeded of e„i in n' independent trials. In the limit of large n' 
and small e this can be approximated as y^: - -2 ln[e„'/«']. 
Of course, if we knew which frequencies to test a priori then 
we could perform < «' independent tests and the significance 
levels could be adjusted accordingly. 

4.2. Accounting for model uncertainty 

The case outlined above is valid only when we know the true 
power spectrum exactly (Pj - V j). In reality all we have is an 
estimated model Pj (which will differ from the true spectrum) 
and its uncertainty. This extra uncertainty alters the probability 
distribution. The ratio yj - 2Ij/Pj is really the ratio of two ran- 
dom variables; the PDF of this would allow us to calculate the 
probability of observing a given value of yj taking full account 
of the uncertainty in the model fitting. As discussed above 2Ij 
will follow a rescaled;^'^ distribution about the true spectrum. 



1 



-x/2Pj 



(17) 



In the case of the LS fitting discussed in section|31the model 
Pj has a log-normal distribution. The probabiUty distribution of 
the power in the fitted model at frequency fj is therefore: 



ppiy) 



Sjyy/2^ 



exp< 



(\n[y]-Mjf 
2S^ 



(18) 



where Mj - ln[!Pj] and S j is the uncertainty on the logarithm 
of the model (equation I13>. The periodogram ordinate Ij and 
the powe Pj in the best fitting model are only strictly inde- 
pendent if the frequency of interest, fj, is excluded from the 
LS fit. Otherwise the fitted model, and hence the model power 
at this frequency, Pj, is influenced by Ij. Although the effect is 
small (~ 1 In) it does have a substantial impact on the tail of the 
PDF. Thus, to obtain formally independent variables we must 
calculate P j (and its error) using the LS method after ignoring 
Ij from the fit. In other words, we must ignore the candidate 
frequency when fitting the continuum model and then compare 
the measured periodogram ordinate at this frequency with the 
model derived from fitting all the other (independent) frequen- 
cies. 

The PDF of the ratio yj can be obtained using the standard 
formula for the PDF of the ratio of two independent variables: 



Py,(z)= r 



\y\P2i,{zy)P'p,(y)dy 



(19) 



The periodogram is always positive so we can integrate over 
positive values only. 



1 



2Sj'Pj^l2^ 



f 

Jo 



exp<^ — 



(ln[ 7] - MY 

2S] 



-^U.(20) 
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The dummy variable w - y/Pj (and dy - Pjdw) can be used to 
simplify the above equation: 



Pyjiz) = 




It is worth noting that this formula contains no dependence on 
the actual value of (and hence Mj). This should not be sur- 
prising because we are dealing with the distribution of the ratio 
of the data to the model, not the absolute value of the data, 
and so the absolute value of the model is not relevant. The 
important parameter is S j, describing the uncertainty on the 
model, and this can be calculated from the theory of linear re- 
gression (section IT^ . Unlike the ideal case discussed above 
(section im the PDF is frequency-dependent, this is reflected 
in the changes of S j with frequency. As S j — > this formula 
reduces to the equivalent for the ideal case discussed above. 

For a given frequency fj the integral in equation|^can be 
evaluated numerically to give the PDF for fj. Figure Q com- 
pares the prediction of equation]^ with a Monte Carlo distri- 
bution at two diff'erent frequencies. Also shown for comparison 
is the X2 PDF which represents the distribution in the absence 
of uncertainties on the model (i.e. 5y — > 0). At small y (i.e. 
low significance peaks) the two distributions agree, whereas for 
large y (i.e. high significance peaks) there is a substantial dif- 
ference between the PDFs including and excluding the uncer- 
tainty on the model. This means that, while the effect of includ- 
ing the uncertainty on the model is negligible for low signifi- 
cance peaks, the significance of high significance peaks may be 
substantially overestimated if this additional uncertainty is not 
taken into account. 

The probability of obtaining a value of jj higher than 
can be computed by integrating this PDF: 

Pr{r;>r.)= Pyj(z)dz=ei (22) 

This can be evaluated numerically to find y^ for a given ei. 
Equivalently, we can find the value of y^ at the corresponding 
1 - ei significance level: 

Pr{f ; < r.) = r PrM^dz = 1 - ei (23) 
Jo 

The calculation of y^ depends only on Py^iz), from equation]^ 
which in turn depends only on S j, from equation[2l and this is 
calculated using the the abscissae (frequencies fj) with no de- 
pendence on the ordinates (periodogram powers 7^). The criti- 
cal value 7e can be evaluated using only the frequencies of the 
periodogram. 

Finally, we need to correct for the number of frequencies 
examined. The probability that a peak will be seen given that 
n' frequencies were examined is e„' - 1 - (1 - ei)" . One can 
find the global (1 - e„')100 per cent confidence level by finding 
the value y^ that satisfies: 

Py^(z)dz = 1 - (1 - e„0'^"' * e«'/n' (24) 

where n' is again the number of frequencies examined. 




Fig. 7. Monte Carlo demonstration of the PDF of the ratio yj. 
The histograms represent the distribution of the ratio yj (mea- 
sured at the j = 10 and j = 100 Fourier frequencies) from 
10^ Monte Carlo realisations of K - 256 time series (with an 
a - 2,N - 1 spectrum). For each of the 10^ simulated time 
series the periodogram was fitted using the LS method, after 
ignoring frequency j, and the ratio yj - lljiPj was measured. 
The solid curve marks the predicted PDF using eauationOTIand 
the dotted line marks the PDF of a.x\ distribution (equation|5}. 
The difference between the two model PDFs is due to the un- 
certainty on the model fit. (Compare with Fig. 5 of Israel & 
Stella il996» . 

As an illustration of the effect of the model uncertainty, 
consider a peak in the j - IQ frequency bin of a n = 128 peri- 
odogram. Neglecting the effect of model uncertainty the nom- 
inal ei - lO"'* threshold is y^ - 18.42 (using equation Il5> . 
But, after including the uncertainty in the model, the proba- 
bility of this level being exceeded is really 3.6 x lO"** (using 
equation l22l . For n - 128 trials this corresponds to global sig- 
nificances of 98.7 per cent confidence (ignoring the model un- 
certainty) and 95.4 per cent confidence (including the model 
uncertainty). The first of these might be called a significant de- 
tection, but once the model uncertainty is taken into account the 
detection is no longer very significant. The difference is even 
more profound for higher significances. 

5. Verification of the method 

The procedures discussed above were tested using Monte Carlo 
simulations. The simulations measured the type I error rate, or 
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Fig. 8. Monte Carlo study of the performance of tests for sig- 
nificant periodogram peaks. For each panel 10^ random time 
series of length K = 256 were generated (with an a = 2, = 1 
spectrum). For each of the 127 frequencies tested {- K/2 - 1 
since the Nyquist frequency was ignored) the fraction of Monte 
Carlo simulations with a peak exceeding the nominal 99.9 and 
99.99 per cent confidence levels was recorded. The upper panel 
used the confidence levels from equation The lower panel 
used the confidence levels computed with equation HI] which 
accounts for the uncertainty in the model. The error bars were 
computed using -^^(1 - p)/N where is the number of simu- 
lations. 



the rate of 'false positive' (spurious) detections of periodic sig- 
nals. For this experiment many artificial time series were gen- 
erated based on a power law spectrum (i.e. the null hypothe- 
sis). For each simulation the threshold, coiTesponding to a 
1 - e confidence level, was calculated and the number of peri- 
odogram ordinates that exceeded this value were recorded. The 
rate measured from the Monte Carlo simulations should be the 
same as the false alarm probability e, often called the 'size of 
the test,' which is the expected rate of type I errors. If the ob- 
served rate of false detections exceeds the nominal size of the 
test then one should expect an excess of spurious detections 



(detections may not be reliable). If the observed rate falls be- 
low the nominal test size then the test is conservative (it gives 
even fewer spurious detections than expected). 

The Monte Carlo rate was derived from 10^ random time 
series of length K - 256 (generated with a or = 2,A^ = 1 
spectrum). For each series the periodogram was computed and 
fitted using the LS method. In the first run, the effects of the un- 
certainty on the model were ignored and the thresholds were 
computed for ei = 10"^' and 10""^ using equation[2] These cor- 
responds to 99.9 and 99.99 per cent confidence levels in a sin- 
gle trial test. For each frequency the fraction of simulations that 
show peaks larger than the threshold was recorded. As shown 
in Fig.|8](upper panel) the observed rate of type I errors in the 
simulated data was far in excess of the nominal size of the test. 
Thus the actual rate of spurious detections was higher than the 
nominal test size, and greatly so at low frequencies where the 
model is more uncertain. The situation is worse at high signif- 
icances (small test sizes) because the tail of the PDF diverges 
from the expectation (Fig.Q. This means that significances cal- 
culated by equation[21will be overestimated. 

The situation is much better when the threshold was 
computed (again for ei = 10"^, 10""^) using equation!^ Again 
10^ random time series were generated with K = 256 using the 
same spectrum. For each periodogram the model power Pj, and 
its error parameters S y, were computed by ignoring frequency 
fj and then fitting with the LS method. The ratio yy was com- 
pared to the critical threshold (computed from equation l22l 
at each frequency. The fraction of Monte Carlo simulations that 
showed 'significant' peaks was very close to the expected level 
(the nominal size) and independent of frequency. The excep- 
tions are the » 5 lowest frequencies. Here the model is least 
constrained and the assumption implicit in equation[n that the 
distribution of log[!Py] is normal, becomes inaccurate. For j>5 
the confidence levels predicted by equation|22]gave the correct 
rate of type I errors. 

Figure |5] shows a specific example, namely the same data 
as in Fig. [2 with the LS power law model. Also shown are 
the ('global') n'-trial confidence limits computed as discussed 
in section l4~2l Clearly none of the peaks in the periodogram 
exceeds the 95 per cent limit, as expected for red noise. 

6. Caveats and comparison with otiier metliods 

6.1. Underlying assumptions 

In order for the test to give reliable significance limits the un- 
derlying noise spectrum must be a power law. Clearly if the 
broad-band noise spectrum does not resemble a power law the 
results of the LS fitting will not be valid. The general solution 
to this problem is to replace the LS fitting procedure with the 
exact maximum likelihood (ML) procedure for fitting the 
distributed periodogram. The appendix describes this method. 

The test was intended to be used for assessing the signif- 
icance of peaks in the periodograms of X-ray observations of 
Seyfert galaxies, which tend to be rather short K ~ 10^ and also 
show significant variance due to measurement errors (Poisson 
noise). These measurement uncertainties produce a flat com- 
ponent that is added to the source power spectrum in the peri- 
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Fig. 9. Same periodogram as in Fig.^ Plotted are the de-biased 
LS estimate of the power law spectral model (solid curve) and 
the 95 and 99 per cent upper limits (dotted curves) on the ex- 
pected power (global significance levels for n - 127 indepen- 
dent frequencies). 



odogram. The effect will cause the observed spectrum to flatten 
at high frequencies as the power in the red noise spectrum of 
the source becomes comparable to the power in the flat Poisson 
noise spectrum. Using the normalisation given in equation^the 
expected Poisson noise level is = 2{{x} + B)/{x}^ where (x) 
is the mean count rate and B is the mean background rate^. It 
is also now known that at low frequencies the power spectra of 
Seyferts break from a single power law (e.g. Uttley et al. 2002; 
Markowitz et al. 2003). These deviations from a single power 
law should be accounted for in modelling the spectrum. The 
simplest solution is to divide the periodogram into frequency 
ranges within which the power spectrum is approximately a 
single power law. The period detection test can then be used as 
described above. The crucial point is that as long as the peri- 
odogram can be fitted reasonably well with a power law over 
a frequency range of interest (as judged using the KS test) the 
test will be valid. Alternatively one may fit a model of a power 
law plus constant (to account for the flattening) using the ML 
method discussed in the Appendix. 

Furthermore, the test, which is based on the discrete Fourier 
transform, is most sensitive to sinusoidal periodicities. Non- 
sinusoidal variations will have their power spread over sev- 
eral frequencies which will lessen the detection significance in 
any one given frequency. Other methods such as epoch folding 
(Leahy et al.[T983 ), whereby one bins the time series into phase 
bins at a test period, can be more sensitive to such variations. 
At the correct period the periodic variations will sum while any 
background noise will cancel out, thus revealing the profile of 
the periodic pulsations. However, the background noise will 
only cancel out if it is temporally independent, i.e. white noise. 
Again, the presence of any underlying red noise variations may 



produce unreliable results (Benlloch et al. l200H if not correctly 
accounted. 

6.2. Comparison with other methods 
6.2.1 . Lomb-Scargle periodogram 

The method described above is based around the standard 
(Fourier) periodogram and therefore requires uniformly sam- 
pled time series. This ensures the asymptotic independence of 
the periodogram ordinates. If the time series is non-uniformly 
sampled one may use other periodogram estimators such as the 



^ For the case of Gaussian errors tr^. on each flux measurement the 
expected Gaussian noise level is = 2Ar(cr'>/(x>^. This assumes 
the time series comprises contiguous bins of length Ar. 



Lomb-Scargle periodogram (Lomb J^976. Scargle [19821 Press 
& Rvbicki ri989> . However, the behaviour of these will not be 
identical to that discussed above. The above procedure should 
not be used on non-uniformly sampled time series (nor should 
the method discussed in section 13.8 of Press et al. 1 19961 be 
used in the presence of non-white noise). Zhou & Sornette 
(2001) discuss the results of Monte Carlo tests on the distri- 
bution of peaks in the Lomb-Scargle periodogram for various 
types of processes. 

6.2.2. Oversampled periodogram 

Oversampling the periodogram, i.e. calculating periodogram 
ordinates at frequencies between the normal Fourier frequen- 
cies, is sometimes done in order to increase the sensitivity to 
weak signals that lie at frequencies between two Fourier fre- 
quencies. A periodic variation with a frequency nearly equidis- 
tant between two adjacent Fourier frequencies, e.g. between fj 
and fj+i, will have its power spread (almost entkely) between 
these two frequencies, thus reducing the significance in any one 
frequency bin. The reduction in power per bin can be as much 
as ^ ^ 0.41. In these situations oversampling the pe- 

riodogram by including additional frequencies between fj and 
fj+i can increase sensitivity to the periodicity. However, it must 
also be noted that by oversampling the periodogram one is test- 
ing more than K/2 frequencies, allowing many more opportu- 
nities to find spurious peaks. The number of trials increases 
above the usual n - Kjl case if the periodogram is oversam- 
pled, although the effective number of independent trials does 
not scale linearly with the oversampling factor because there 
is a fixed number in) of strictly independent frequencies (the 
Fourier frequencies). 

Oversampling the periodogram and assuming ^ n trials 
will tend to overestimate the significance of peaks in the peri- 
odogram. This was demonstrated by Monte Carlo simulations 
(see Fig. not . For this demonstration 10"^ white noise time series 
(spectrum: a - Q,N - 1) were generated with length K - 256 
and for each one the periodogram was calculated using both 
the standard Fourier frequencies and also oversampling the fre- 
quency resolution by a factor 8. The peak power in the peri- 
odogram of each of the simulated time series was recorded. 
The distribution of peak powers is shown in Fig.QSIfor both the 
standard and the oversampled periodograms. There is an obvi- 
ous tendency for the oversampled periodogram to show slightly 
higher peak values, as might be expected based on the above 
arguments. These peak powers were translated to 'global' sig- 
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nificances (assuming n independent trials) and the distribution 
of significances is shown in the lower panel of Fig. [TO] The 
distribution is flat for the Fourier sampled periodogram: the 
significance of the peaks is exactly as expected. The distribu- 
tion derived from the oversampled periodogram clearly shows 
a substantial excess of high significance peaks. This means that 
the oversampled periodogram is likely to produce many more 
spurious peaks than expected if one assumes only n indepen- 
dent trials were made. A Monte Carlo estimate of the global 
significance would be required to calibrate the significance of 
peaks in oversampled periodogram (and find a more realistic 
effective number of trials). 

An alternative to oversampling the periodogram is to per- 
form a sliding two-bin search for periodogram peaks using 
the standard Fourier frequencies only, as discussed by van der 
Klis y_989 section 6.4). This increases the sensitivity to pe- 
riods whose frequencies fall between the Fourier frequencies. 
However, one will still need to perform Monte Carlo simula- 
tions to assess the global significance (by measuring the rate of 
type I errors in the simulations). 

6.2.3. Monte Carlo testing 

An alternative test for periodic variations in red noise is to es- 
timate the likelihood of observing a given peak using Monte 
Carlo simulations of red noise processes (Benlloch et al. 2001 
Halpern, Leighly & Marshall l2003l l. However, the method of 
Benlloch et al. (i2001i does not account for uncertainties in the 
best-fitting model, which can seriously effect the apparent sig- 
nificances of strong peaks (section RT^ . Monte Carlo simula- 
tions are only as good as the model they assume! One solution 
would be to map the (multi-dimensional) distribution of the 
model parameters and for each simulation draw the model pa- 
rameters at random from this distribution. This would thereby 
account for the uncertainty in the model parameters. (Protassov 
et al. l2002l use a similar approach to calculate posterior predic- 
tive p-values.) 

7. Discussion 

A simple procedure is presented for assessing the significance 
of peaks in a periodogram when the underlying continuum 
noise has a power law spectrum. 

7.1. Recipe 

The following is one possible recipe for periodogram analysis. 

- Calculate the periodogram of the data I(fj) (section^ and 
convert to log-space, i.e. log[fj] and log[/(/j)]. 

- Ignore the Nyquist frequency and frequencies above which 
Poisson noise is significant, leaving «' frequencies. 

- Estimate the power law parameters by fitting a linear func- 
tion to the log-periodogram using the LS method (sec- 
tion|3}. 

- Test the goodness of the fit by comparing the distribution 
of data/model residuals with the expectation using a KS 
test (section|3}. 




Peak power 
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Fig. 10. Comparison between periodograms sampled at the 
K/2 Fourier frequencies and oversampled by a factor 8. The 
upper panel shows the distribution of peaks in the periodogram 
from 10"* random white noise time series of length K = 256 (us- 
ing an a = 0, = 1 spectrum). For each simulation the value 
of the peak of the periodogram was recorded for the standard 
(solid histogram) and oversampled (dotted histogram) case. 
The oversampled periodograms clearly show slightly larger 
peak values. The lower panel shows the same data except the 
peak powers have been converted to 'global' significances us- 
ing eauation[T6l 



- Calculate the 'global' («'-trial) threshold ratio for a (1 - e) 
significance detection using = -2 ln[l -(1 -e)'^" ] where 
n' is number of periodogram points used in the fit and e is 
the desired 'false alarm probability' (section lTTl . 

- Multiply the best-fitting continuum model by (or add 
its logarithm in log-space). 

This will provide an estimate of the power spectral slope 
and normalisation and a first indication of the presence of sig- 
nificant periodicities. However, as discussed in section 14^21 the 
uncertainty inherent in the model fitting means that the sig- 
nificances of strong peaks may be substantially overestimated. 
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This procedure should be quite accurate for low significance 
peaks however, meaning that the above procedure is valid for 
rejecting low significance peaks (ej > 0.01). The inaccuracy 
of the method at high powers means that the procedure should 
only be used to reject low significance peaks, not detect high 
significance peaks. In order to obtain a more rigorous estimate 
of the significance of stronger peaks one needs to account for 
the additional uncertainty in the model. 

For a given frequency of interest, fj, one must: 

- Ignore fj from the fit (in order that the data and model are 
independent at fj). 

- Re-fit the periodogram using the LS method (section|3}. 

- Calculate the power in the model at this frequency, Pj. 

- Measure the ratio yj - lljjPj. 

- Compute the uncertainties on the model parameters (sec- 
tion l3.2l i and hence the uncertainty on the model continuum 
S j ("section lXSt . 

- Calculate the probability of observing a value of jj this 
high by numerically integrating equationl22l 

The above procedure is essentially the application of two 
standard tools for time series analysis. The first is estimat- 
ing the power law spectral properties by LS fitting of the 
log-periodogram. This has been discussed by several authors 
(e.g. Geweke & Poi-tei--Hudak .l983. Fougere [T985l Pilgram & 
KaplanQWSJ. The second tool is applying the known ;^^2 prop- 
erties of the periodogram to estimate confidence levels, as is 
standard for white noise spectra (Priestlev .1981. van der Klis 
I1989> . The additional calculation to include the uncertainty in 
the model follows standard statistical procedures. The overall 
approach is similar to that discussed by Israel & Stella U996 ) 
but is tailored to power law spectra. 

In order for the method to produce reliable results the un- 
derlying power spectrum has to be a power law. However, the 
test will work well even for data that show deviations from a 
power law (such as intrinsic low frequency flattening or high 
frequency flattening due to Poisson noise) provided the peri- 
odogram is divided into frequency intervals within which a sin- 
gle power law provides a good description of the data (deter- 
mined with a KS test). In this case the power spectrum over the 
restricted frequency ranges is indistinguishable from a power 
law and, when applied to these limited frequency ranges, the 
method will function as expected. 

7.2. Application to real data 

As a demonstration of the method we present a re-analysis 
of the X-ray observations of two Seyfert galaxies. The first is 
the long ASCA observation of IRAS 18325-5926 (Iwasawa et 
al. [1998) and the second is the XMM-Newton GTO observa- 
tion of Mrk 766 (Boiler et al. l2001t . For IRAS 18325-5926 a 
background-subtracted 0.5 - 10 keV SISO light curve was ex- 
tracted in 100-s bins, using a 4 arcmin radius source region, 
and rebinned onto an evenly spaced grid at the spacecraft or- 
bital period (5760-s). For Mrk 766 a background-subtracted, 
exposure-corrected 0.2 - 2 keV light curve was extracted from 
the EPIC pn using a 38 arcsec radius source region. The re- 
binned IRAS 18325-5926 hght curve contained 84 regularly 
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Fig. 11. Periodogram of IRAS 18325-5926 from ASCA. 
Plotted are the de-biased LS estimate of the power law spec- 
tral model (solid curve) and the 95 and 99 per cent upper limits 
(dotted curves) on the expected power (global significance lev- 
els for n = 19 independent frequencies). Also shown is the 
expected level of the Poisson noise power. 

spaced bins and the Mrk 766 light curve contained 330 con- 
tiguous 100-s bins. 

The periodogram of each light curve was calculated and 
the expected Poisson noise level determined. Only the low- 
est frequency periodogram points were examined, to minimise 
the effect of the Poisson noise level. For IRAS 18325-5926 
only the 19 lowest frequency points were used, while for Mrk 
766 only 97 points were used. These were fitted with a linear 
function in log-space, giving slopes of a ^ 0.9 and a ss 1.8 
for IRAS 18325-5926 and Mrk 766, respectively. The 95 and 
99 per cent confidence limits on the periodogram were com- 
puted, accounting for the number of frequencies examined in 
each case (see figures fTTI and I12>. In neither object did a sin- 
gle periodogram point exceed the 95 per cent limit, meaning 
that there is no strong evidence to suggest a periodic com- 
ponent to their variability, contrary to the original claims of 
Iwasawa et al. a998, f? * 1.7 x 10"^ Hz) and Boiler et al. 
QOOl , fp ^ 2.4 X 10"^ Hz). Benlloch et al. (200 Ij drew similar 
conclusions based on Monte Carlo simulations. 

7.3. Conclusions 

The simplicity of the proposed test means that it can be used as 
a quick, first check against spurious periodogram peaks, with- 
out the need for extensive Monte Carlo simulations. Examples 
where this test might be useful include not only X-ray obser- 
vations of Seyferts but also testing for periodic components in 
monitoring of blazars (e.g. Kranich et al. 119991 Hayashida et 
al. 19981, the Galactic centre (Genzel et al. I2003> and other 
astrophysical sources that show red noise variations. 

The basic idea behind the method is to de-redden the pe- 
riodogram by dividing out the best-fitting power law and then 
using the known distribution of the periodogram ordinates to 
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estimate the likelihood of observing a given peak if the null 
hypothesis (power law spectrum with no periodicity) is true. 
However, this should not be treated as a 'black box' solution to 
the problem of detecting periodicities in red noise data. There 
is no substitute for a thorough understanding of the nuances of 
power spectral statistics, as illustrated by the PDFs discussed in 
section 14^21 The tail of the PDF is very sensitive to the details 
of the method, and thus one needs to treat any periodogram 
analysis with great care in order to avoid overestimating the 
significance of a peak. 

To date the most significant candidate periodicities found 
in the X-ray variations active galaxies are those in long EUVE 
observations of three nearby Seyferts by Halpern et al. (.2003i) . 
However, their published global significances will be overesti- 
mates for two reasons: (/) no account was made of the uncer- 
tainty on the model and (ii) the periodograms were oversam- 
pled. These two separate issues both act to artificially boost the 
apparent significance of spurious signals, as discussed in sec- 
tions |^] and respectively. Similarly the apparently sig- 
nificant candidate periodicity in NGC 5548 found by Papadakis 
& Lawrence (ll993bt was shown to be substantially less sig- 
nificant when the uncertainty in the modelling was included 
(Tagliaferri et al. 119961 ). This leaves us in the interesting situa- 
tion of there being not one surviving, robustly determined and 
significant (> 99 per cent) periodicity in the X-ray variability 
of an active galaxy. 

Finally we close with a plea. There exist many more AGN 
light curves and periodograms than have been published. There 
is a natural publication bias: those that show the most 'signif- 
icant' features get published. This means the true number of 
trials is much larger than for any one given experiment. In such 
situations one should treat with caution the significance of the 
subset of results that are published (e.g. Scargle 21)001. That is, 
until a result is published that is so significant that it cannot 
be accounted for by publication bias. The detection of peri- 
odic or quasi-periodic variations in galactic nuclei would be a 
major discovery and of great importance to the field. The im- 
portance of the result should be argument enough for very high 
standards of discovery. We would therefore advocate serious 
further investigation of only those candidate periodicities with 
high significances (such as a > 99.9 per cent confidence, or a 
"3cr minimum," criterion) after accounting for all likely sources 
of error 
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Appendix A: Maximum liltelihood (l\1L) fitting 

In this section we briefly elucidate the exact method of max- 
imum likelihood (ML) fitting periodograms. See Anderson, 
Duvall & Jeff'eries (1990 1 and Stella et al. (1997 i for more 
details. The probabiUty density function (PDF) of the peri- 
odogram ordinates at each frequency are given by: 
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Fig. 12. Periodogram of Mrk 766 from XMM-Newton. 

where Pj is the true underlying spectrum (the expectation value 
of the periodogram) at frequency fj. This is valid for frequen- 
cies j = 1, 2, . . . , n - 1 since the periodogram ordinate at the 
j = n (Nyquist) frequency is distributed differently. 

Assuming a model Pj{9k), determined by parameters Ok - 
{9\,92, . . . , 9m], we can write the joint probability density of 
observing the n-1 periodogram ordinates: 

X=fj ;,(/,). fj^e-V^^ (A.2) 

./=i i=i 

If the data Ij have already been observed this represents the 
likelihood function. Maximising the likelihood £. is equivalent 
to minimising S = -2 ln[X] and we can rewrite this as: 

5 =2U|ln[!PJ + ^| (A.3) 

Finding the model parameters Ok that minimise S will yield the 
maximum likelihood parameter values. 

One can then use standard tools of maximum likelihood 
analysis, such as the likelihood ratio test (LRT) to test for addi- 
tional free parameters in the model: 

R^-2\n[£il£2\^S,-S2 (A.4) 

where Xi represents the likelihood for the simpler (more parsi- 
monious) model and £.2 represents the likelihood for the model 
with the additional free parameters. The second, more complex 
model, always contains within it the simpler model as a sub- 
set (the models are nested) and therefore the UkeUhood for the 
more complex model is always equal to or greater than that for 
the simpler model. With certain restrictions (see Protassov et 
al. l2002l l R is distributed as a.xl variable where v is the number 
of additional free parameters. This test could be used, for ex- 
ample, to compare nested models such as broken and unbroken 
power laws. 

Alternatively one can compare different models using the 
Akaike Information Criterion (A/C; Akaike . 19731 : 

Aid = -2 ln[i;;] + 2ki = 5 , + 2ki (A.5) 
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where X/ is the likelihood and fc, is the number of free pa- 
rameters for model /. The second term in the sum is a penalty 
for including more free parameters. The model that minimises 
the AIC is considered to be the best (the models need not be 
nested). 

One can also use AS = -2A \n[X.] to place confidence lim- 
its on the model parameters in a fashion exactly analogous 
to mapping confidence contours using Ax^ (section 15.6 of 
Press et al. fT996). Under fairly general conditions (see Cash 
I1979> . e.g. the ln[i]]-surface is approximately shaped like a 
multi-dimensional paraboloid, AS distributed as xl where v 
is the number of parameters of interest (e.g. v = 1 for the 
one-dimensional confidence region on an individual parame- 
ter). One can use standard tables of xl values to place confi- 
dence limits (e.g. AS = 2.71 corresponds to 90 per cent confi- 
dence limits on one parameter). 

For the purposes of period searching one may fit a suitable 
M-parameter continuum model (representing the null hypothe- 
sis, i.e., no periodic signal) using the ML method and define the 
M-dimensional distribution of its parameters (using AS). One 
can then use this M-dimensional distribution of the model pa- 
rameters to randomly draw models for Monte Carlo simulation. 
This procedure will thereby account for the likely distribution 
of model parameters (which gives rise to uncertainties in the 
estimated continuum level). 
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